#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os
import sys
import math
import random

caledDAP=[]
lc=[]

def getData(filepath):
        f=open(filepath,"r")
        pg=f.read()
        f.close()

        data=[]
        s=pg.split("\n")
        for item in s:
                tmp=item.split("\t")
                dd=["%s"%tmp[0]]
                dd.append([float(j) for j in tmp[1:]])
                data.append(dd)
        return data

def calCentralPoint(cp):
	'''cp like this:
	[[xx,[point1]],,[xx,[point2]]...[xx,[pointN]]]
	'''
	length=len(cp)
	if length==0:
		return [0,0]
	cptmp=[]
	for item in cp:
		cptmp.append(item[1])
	tmp=[[row[col] for row in cptmp] for col in range(len(cptmp[0]))]
	center=[]
	for index in range(len(tmp)):
		totals=sum(item for item in tmp[index])
		result=totals/length
		center.append(result)
	return center


def calDistance(x,y,p2p=1):
	if p2p==1:
		return math.sqrt(sum(pow(x[1][index]-y[1][index],2) for index in range(len(x[1]))))
	else:
		return math.sqrt(sum(pow(x[1][index]-y[index],2) for index in range(len(y))))

def kmeans(data,kvalue=2):
	clusters=[[] for j in range(kvalue)]
	centers=[]
	#sampling kvalue numbers in the range of 0~length of data
	randcenterlist=random.sample(xrange(len(data)),kvalue)
	#using the sampled numbers as index to randomly form k clusters
	for index in range(kvalue):
		clusters[index].append(list(data[randcenterlist[index]]))
		centers.append(list(data[randcenterlist[index]][1]))
	#remove points which was put in the clusters from data
	rmdata=[data[randcenterlist[j]] for j in range(kvalue)]
	for j in range(kvalue):
		data.remove(rmdata[j])
	#divide data into k clusters according to the randomly sampled data
	while len(data)!=0:
		minDistance=100000000
		minIndex=0
		#calculate distance between data[0] and the center of every cluster		
		for index in range(kvalue):
			d=calDistance(data[0],centers[index],0)
			#get the minimum distance,
			if d<minDistance:
				minDistance=d
				minIndex=index
		#put the item into its cluster
		clusters[minIndex].append(list(data[0]))
		#recalculating the center of the cluster
		centers[minIndex]=list(calCentralPoint(clusters[minIndex]))
		#remove the classified item
		del data[0]
	return clusters

def chooseMaxCluster(clusters):
	if len(clusters)==0:
		return -1
	maxClusterLength=0
	maxClusterIndex=-1;
	for index in range(len(clusters)):
		if len(clusters[index])>maxClusterLength:
			maxClusterLength=len(clusters[index])
			maxClusterIndex=index
	return maxClusterIndex

def calClusterRadius(cluster):
	'''cluster should like this:
	[[point1],[point2]...[pointN]]
	'''
	center=calCentralPoint(cluster)
	clusterRadius=0
	for item in cluster:
		d=calDistance(item,center,0)
		if d>clusterRadius:
			clusterRadius=d
	return clusterRadius

def stageone(filepath):
	data=getData(filepath)
	#nc is an empirical value
	nc=max(2,len(data)/100)
	kClusters=[]
	#print len(data)
	while len(kClusters)<=nc:
		tmp=kmeans(data)
		for item in tmp:
			kClusters.append(item)
		if len(kClusters)==nc:
			break
		rmindex=chooseMaxCluster(kClusters)
		#data changes to the max clusters,and it will be divided
		data = kClusters[rmindex]
		del kClusters[rmindex]
	odcList=[]
	for item in kClusters:
		rc=calClusterRadius(item)
		c=len(item)
		odcList.append(rc/c)
	#sorting k clusters using odc sequence
	lsttobesort=[(odcList[j],j) for j in range(len(odcList))]
	lsttobesort.sort()
	#descending order
	lsttobesort.reverse()
	resultClusters=[]
	for item in lsttobesort:
		resultClusters.append(kClusters[item[1]])

	return resultClusters


def neighbors(p,collection,k):
	neighborList=[(calDistance(p,item,1),item) for item in collection]
	neighborList.sort()
	if k<len(neighborList):
		return [item[1] for item in neighborList[0:k]]
	return [item[1] for item in neighborList]



def calDAP(p,S,k):
	global caledDAP
	global lc
	if len(S)<=1:
		return 0
	'''
	for item in caledDAP:
		if item[0]==p and S==lc[item[1]]:
			return item[2]
	'''
	ret=0
	index=0
	neighbors=[calDistance(p,item2,1) for item2 in S if p!=item2]
	neighbors.sort()
	if k<=len(neighbors):
		ret=sum(neighbors[0:k])/k
	else:ret=sum(neighbors)/len(neighbors)
	'''
	for j in range(len(lc)):
		if S==lc[j]:
			index=j
	tmp=(p,index,ret)
	caledDAP.append(tmp)
	'''
	return ret


def maxDistance(point,neighborList):
	if len(neighborList)==0:
		return 100000000
	return max(calDistance(point,item,1) for item in neighborList)

def topoutliers(S,n,k):
	if len(S)<=1:
		return S
	outliers=[]
	itemdap=0
	for item1 in S:
		item1neighbors=[calDistance(item1,item2,1) for item2 in S if item1!=item2]
		item1neighbors.sort()
		if k<=len(item1neighbors):
			itemdap=sum(item1neighbors[0:k])/k
		else:itemdap=sum(item1neighbors)/len(item1neighbors)
		outliers.append((itemdap,item1))
	outliers.sort()
	outliers.reverse()
	if n<len(outliers):
		return [item[1] for item in outliers[0:n]]
	return [item[1] for item in outliers]

def calomin(S,k):
	if len(S)<=1:
		return 0
	dap=[]
	itemdap=0
	for item1 in S:
		item1neighbors=[calDistance(item1,item2,1) for item2 in S if item1!=item2]
		item1neighbors.sort()
		if k<=len(item1neighbors):
			itemdap=sum(item1neighbors[0:k])/k
		else:itemdap=sum(item1neighbors)/len(item1neighbors)
		dap.append(itemdap)
	return min(dap)

def stagetwo(clusters,n,k):
	global lc
	lc=clusters
	omin=0
	outlier=[]
	jump=False
	goon=False
	for c in lc:
		for p in c:
			nnp=[]
			for ac in lc:
				for q in ac:
					if p==q:
						continue
					if calDistance(p,q,1)<omin-calDAP(q,ac,k):
						print "%s %-19s %s"%("Cut by #2:",p[0],p[1],)
						jump=True
						goon=True
						break
					if len(nnp)<k or calDistance(p,q,1)<maxDistance(p,nnp):
						nnp.append(q)
						nnp=neighbors(p,nnp,k)
						jump=False
						goon=False
					if len(nnp)==k and calDAP(p,nnp,k)<omin:
						print "%s %-19s %s"%("Cut by #1:",p[0],p[1],)
						jump=True
						goon=True
						break
				if jump==True:
					break
			if goon ==True:
				continue
			print "%s %-19s %s"%("Outlier  : ",p[0],p[1],)
			outlier.append((calDAP(p,nnp,k),p))
			outlier.sort()
			outlier.reverse()
			if len(outlier)>n:
				outlier=outlier[0:n]
			omin=min(item[0] for item in outlier)
			#print "omin changed=======================","new omin",omin
	print
	print
	print "结果:   共有%s个离群点:"%(len(outlier),)
	print
	print "%-8s  %-19s      %s"%("离群因子","身份证号","数据",)
	for item in outlier:
		print "%-8.2f  %-19s  %s"%(item[0],item[1][0],item[1][1])
	return



stagetwo(stageone("e:/yibao1.txt"),30,30)
  		
  		
